Simulating A Pendulum's Motion¶
Inspired largely by this video on differential equations by 3Blue1Brown. In this notebook, we model the motion of a pendulum using a differential equation and then simulate it using Python scientific computing tools.
import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
import matplotlib.animation as anim
from IPython.display import HTML
%matplotlib notebook
The Pendulum Equation¶
A Free Body Diagram, or FBD, outlines all of the forces acting on an object, and it better helps understand the system.
Analyzing this FBD, we see that the pendulum is influenced by two forces. First, there's gravity, which always acts downward. However, we understand that the pendulum can only move across the arc defined by the angle $\theta$. So, the only component of the force that influences the pendulum is along the direction of this arc. This, as it turns out, is proportional to $\sin{\theta}$.
So our first force is $F_g \sin{\theta}$, where $F_g = mg$.
The second component is the drag force $F_D$. This is usually modeled as proportional to the velocity of the system. However, velocity can be modeled as a function of the change in the angle (or angular velocity): $v = L \dot{\theta}$. We can define a coefficient $b$ to be the drag coefficient of this system. So the final equation for drag force is $F_D = bL\dot{\theta}$
We also need to understand that the force of gravity acts in the opposite direction of the angle $\theta$, trying to pull our pendulum down closer to the ground. So $F_g$ is actually negative. Likewise, drag force acts in the opposite direction to velocity (moving with velocity causes friction with oncoming air, slowing our pendulum down). So $F_D$ is also negative.
Our final force equation is the follow:
$$ \Sigma F = - mg \sin{\theta} - bL \dot{\theta} $$Newton's Second Law states that the net force on an object is equal to the mass times the acceleration. Acceleration is modeled as proportional to the change in angular velocity $\ddot{\theta}$ by the length $L$. The equation then becomes the following
$$ m L \ddot{\theta} = - mg \sin{\theta} - bL \dot{\theta} $$Rearranging, we get our final equation
$$ \ddot{\theta} = - \frac{g}{L} \sin \theta - \frac{b}{m} \dot{\theta} $$Expand it to a vectorized form (using $\omega = \dot{\theta}$ as angular velocity)
$$ \begin{bmatrix} \dot{\theta} \\ \dot{\omega} \end{bmatrix} = \begin{bmatrix} \omega \\ -\frac{g}{L} sin(\theta) -\frac{b}{m} \omega \end{bmatrix} $$We can now write this equation as a python function which will help our analysis
def pendulum(t, Th, g=9.81, L=1, m=3, b=0.5):
"""
Pendulum differential equation
"""
th, om = Th
omd = - (b/m) * om - (g/L) * np.sin(th)
out = np.array([ om, omd ])
return out
Phase Plot of Equation¶
Differential equations can be better analyzed in phase space. A phase space is an abstract space where the dimensions are the position and velocity components of the system. Each position in the space represents a unique state that the pendulum could be in. The model equation is represented as a vector valued function that shows the direction that the pendulum system will travel next in this space. Given a position in the space, the phase diagram shows the next position, and then the next position after that. So the plot will be a vector plot with the x axis representing the angle of the pendulum and the y axis representing the angular velocity of the pendulum.
thr, omr = np.meshgrid(
np.linspace(-2*np.pi, 2*np.pi, 30),
np.linspace(-2*np.pi, 2*np.pi, 30)
)
thv, omv = pendulum(0.0, [thr, omr])
plt.figure()
plt.quiver(thr, omr, thv, omv, np.hypot(thv, omv))
plt.xlabel('$\\theta$ (radians)')
plt.ylabel('$\omega$ (rads/sec)')
plt.title('Phase Plot of Pendulum Equation')
plt.axis([ -2*np.pi, 2*np.pi, -2*np.pi, 2*np.pi ])
plt.show()
Generate Phase Flows using scipy.integrate¶
# Equation meta parameters
L = 10
g = 9.81
m = 3
b = 0.5
# Time span
t = (0, 20)
max_dt = 0.01
# Initial parameters
Th0s = np.radians([
[ 30, 0 ],
[ 45, 0 ],
[ 60, 0 ],
[ 90, 0 ],
[ 150, 0 ],
[ 190, 0 ],
[ 0, 30 ],
[ 0, 45 ],
[ 0, 60 ],
[ 0, 90 ],
[ 0, 150 ],
[ 0, 180 ],
[ 0, 270 ]
])
# Do all the integrations
Ths = [
spi.solve_ivp(
pendulum,
t, Th0,
max_step=max_dt,
args=(L, g, m, b)).y
for Th0 in Th0s ]
# Plot them all!
plt.figure()
for Th in Ths:
plt.plot(Th[0], Th[1])
plt.axis([-2*np.pi, 2*np.pi, -2*np.pi, 2*np.pi])
plt.xlabel('$\\theta$ (radians)')
plt.ylabel('$\\omega$ (rads/sec)')
plt.title('Flow Plot of Pendulum')
plt.show()
Simulate Multiple Pendulums¶
# Initial parameters
Th0 = np.radians([
# Initial angles
[ 30, 60, 90, 190, 0, 0, 0, 0, 0],
# Initial angular velocities
[ 0, 0, 0, 0, 30, 60, 90, 180, 270]
])
# Equation meta parameters
L = 10
g = 9.81
m = 3
b = 0.5
# Figure parameters
interval = 20
framerate = 30
size = 15
# Turning interactive plot off temporarily
plt.ioff()
# Plot dimensions
dims = 8, 6
ylim = -size, size
xsize = size * dims[0] / dims[1]
xlim = -xsize, xsize
# Create figure
fig = plt.figure(figsize=dims)
ax = plt.axes(xlim=xlim, ylim=ylim)
lines = [
ax.plot([], [], lw=1, marker='o')[0]
for i in range(Th0.shape[1])
]
plt.xticks([])
plt.yticks([])
plt.title('Pendulum Simulation!')
def init():
"""
Initialization function
"""
for line in lines:
line.set_data([], [])
return lines
Th = Th0
def animate(i):
"""
Animation function
"""
# Define our global pendulum variable
global Th
dt = 1/framerate
# Update pendulum state
Th += pendulum(t, Th, g, L, m, b)*dt
# Get x and y components
th = Th[0,:]
x = L*np.sin(th)
y = -L*np.cos(th)
# Draw each pendulum
for i,line in enumerate(lines):
line.set_data([0, x[i]],[0, y[i]])
# Return updated lines
return lines
# Generate animation
frames = interval*framerate
anm = anim.FuncAnimation(
fig, animate,
init_func=init,
interval=interval,
frames=frames,
blit=True)
# Close plot
plt.close(fig)
plt.ion()
# Render animation
HTML(anm.to_jshtml())